
function [exu,Tk,Tkk,Tkw,Ty,Tyk,Tyw,Tkv] = qzs(varargin)
% qzs.m
%
% [exu,Tk,Tkk,Tkw,Ty,Tyk,Tyw,Tkv]
%     = qzs(A1,A2,B1,B2,C1,C2,gam,gamW,gamV,phi,phiW,phiV)
%
% takes a linear rational expectations system of the form
%
%     A1*y(t+1) = B1*y(t) + C1*x(t+1)
%     A2*E(y(t+1)|t) = B2*y(t) + C2*E(x(t+1)|t)
%
%     x(t) = gam + gamW*w(t-1) + gamV*v(t)
%     w(t) = phi + phiW*w(t-1) + phiV*v(t)
%
% and computes a state space solution of the form
%
%     k(t) = Tk + Tkk*k(t-1) + Tkw*w(t-1) + Tkv*v(t)
%     y(t) = Ty + Tyk*k(t) + Tyw*w(t)
%
% Will not handle unit roots. Closely follows appendix for clarity. 
%
% exu is a vector where there the first element is true if a solution
% exists, and the second element is true if the solution is unique.
%
% Options:
% -  'quiet=': logical indicating whether output should be printed to the
%    command window (default: false)
% -  'useoldqz=': keep qz decomposition from previous call, do not
%    re-compute (default: false)
% -  'check=': logical indicating whether to check the existence/uniqueness
%    conditions (default: true, unless existing qz decomposition is being 
%    used, in which case the default is false). 

persistent Q Z N J ns nu R Rs Ru Atilde JbuN Nss Nuu Nsu Jss Juu Jsu

% first parse options
% parseOptions takes inputs: (varargin list, option name,
% and default value if not found)
[quiet,varargin] = parseOption(varargin,'quiet',false) ;
[useoldqz,varargin] = parseOption(varargin,'useoldqz',false) ;
[check,varargin] = parseOption(varargin,'check',not(useoldqz)) ;
if ~useoldqz 
    Q = [] ;
    Z = [] ;
    N = [] ;
    J = [] ;
    ns = [] ;
    nu = [] ; 
    R = [] ;
    Rs = [] ;
    Ru = [] ; 
    Atilde = [] ; 
    JbuN = [] ;
    Nss = [] ; 
    Nuu = [] ;
    Nsu = [] ; 
    Jss = [] ; 
    Juu = [] ;
    Jsu = [] ;
end

% then pull remaining input arguments
% first block of equations holds exactly
A1 = varargin{1} ;
B1 = varargin{3} ;
C1 = varargin{5} ;
% second block of equations holds in expectation
A2 = varargin{2} ;
B2 = varargin{4} ;
C2 = varargin{6} ;
it = 6 ;

nf = size(A1,1) ;
ng = size(A2,1) ;
ny = ng+nf ;
nx = size(C1,2) ;

if numel(varargin)>it
    gam = varargin{it+1} ;
    gamW = varargin{it+2} ;
    gamV = varargin{it+3} ;
    phi = varargin{it+4} ;
    phiW = varargin{it+5} ;
    phiV = varargin{it+6} ;
else
    % use dummy driving process which assumes x(t) is i.i.d. and
    % conformable with C
    gam = zeros(nx,1) ;
    gamW = zeros(nx,1) ;
    gamV = eye(nx,nx) ;
    phi = 0 ;
    phiW = 0 ;
    phiV = zeros(1,nx) ;
end

assert(ny == size(A1,2),'Dimension mismatch.') ;
assert(isempty(A2) || ny == size(A2,2),'Dimension mismatch.') ;
assert(size(B1,1)+size(B2,1) == ny,'[A1;A2] and [B1;B2] must be the same size.') ;
assert(size(B1,2) == ny,'Dimension mismatch.') ;
assert(isempty(B2) || size(B2,2) == ny,'Dimension mismatch.') ;
assert(size(C1,1)+size(C2,1) == ny,'Dimension mismatch.') ;
assert(nx == size(gam,1),'Dimension mismatch.') ;
assert(size(gam,2) == 1,'Dimension mismatch.') ;
assert(size(gamW,2) == size(phi,1),'Dimension mismatch.') ;
assert(size(gam,1) == size(gamW,1),'Dimension mismatch.') ;
assert(size(gam,1) == size(gamV,1),'Dimension mismatch.') ;
assert(size(phi,1) == size(phiW,1),'Dimension mismatch.') ;
assert(size(phi,1) == size(phiV,1),'Dimension mismatch.') ;
assert(size(phiW,1) == size(phiW,2),'Dimension mismatch.') ;

if ~useoldqz
    % Use QZ to create system in canonical variables for decoupling.
    %
    % For reference:
    %     Q*A*Z = N, Q*B*Z = J.
    %     A = Q'*N*Z'
    %     B = Q'*J*Z'
    [N, J, Q, Z] = qz([A1;A2],[B1;B2]) ;
    eigv = ordeig(N,J) ;
    stab = abs(eigv) > 1 ;
    [N, J, Q, Z] = ordqz(N,J,Q,Z,stab) ; % sorts stable roots to upper left block
    
    % count number of stable/unstable roots
    ns = sum(stab) ;
    nu = ny - ns ;
    
    % convert qz output to form used in paper
    Q = Q' ; % inversion because Q is unitary
    Z = Z' ; % inversion because Z is unitary
    R = Z' ;
    Rs = R(:,1:ns) ;
    Ru = R(:,ns+1:ny) ;
    
    Atilde = A1*Rs ;
    
    Nss = N(1:ns,1:ns) ;
    Nuu = N(ns+1:ny,ns+1:ny) ;
    Nsu = N(1:ns,ns+1:ny) ;
    Jss = J(1:ns,1:ns) ;
    Juu = J(ns+1:ny,ns+1:ny) ;
    Jsu = J(1:ns,ns+1:ny) ;

    JbuN = Juu\Nuu ;
end

nw = size(phiW,1) ;

% partition canonical decoupled system
PSI = Q'*[C1;C2] ;
PSIs = PSI(1:ns,:) ;
PSIu = PSI(ns+1:ny,:) ;

% calculate solution for N2tilde in the equation
%     N2tilde - (Jss\Nuu)*N2tilde*phiW = (Juu\Nuu)*(Juu\PSIu)*gamW*(eye(nw)-phiW)\phiW
% using findvss()
%     H - A*H*C = B
% for mapping structural shocks v(t) to the persistent exogenous process w(t)
if max(max(abs(PSIu*gamW)))>eps
    if max(max(abs(phiW)))>eps
        N2tilde = findvss( JbuN, (JbuN)*(Juu\PSIu)*gamW*((eye(nw)-phiW)\phiW), phiW ) ;
    else
        N2tilde = zeros(nu,nw) ;
    end
    Ntilde = (eye(nu)-JbuN)\(JbuN)*(Juu\PSIu)*gamW/(eye(nw)-phiW) - N2tilde ;
else
    Ntilde = zeros(nu,nw) ;
end

% find gUW using H = findvss(A,B,C) where
%     H - A*H*C = B
%     gUW - Juu\Nuu*gUW*phiW = Juu\PSIu*gamW
% for mapping persistent exogenous variables w(t) to canonical
% unstable variables u(t)
if max(max(abs(phiW)))<eps
    gUW = -Juu\PSIu*gamW ;
else
    gUW = -findvss( JbuN, Juu\PSIu*gamW ,phiW ) ;
end
gU = -(eye(nu)-JbuN)\(Juu\PSIu)*gam -  Ntilde*phi ;
gSW = Nss\( Jsu*gUW + PSIs*gamW - Nsu*gUW*phiW ) ;
gS = Nss\( Jsu*gU + PSIs*gam - Nsu*(gU + gUW*phi) ) ;

% existence
if check
    Ss = NaN(nf,nx*(nu+1)) ;
    Ss(:,1:nx) = C1 ;
    for ii = 1:nu
        Ss(:,1+nx*ii:nx*(ii+1)) = A1*Ru*(JbuN^(ii-1))*(Juu\PSIu) ;
    end
    [ex,q,tol] = span(Atilde,Ss) ;
    
    if ~quiet
        if ex
            fprintf(1,'Solution exists; spanning condition met. \n') ;
        else
            fprintf(1,'Spanning condition not met by %g with tolerance %g. \nEither a solution does not exist or some of the equations which \nare specified to hold exactly only hold in expectation.\n',q,tol) ;
        end
    end
    
    % uniqueness
    if rank(Atilde) == ns
        ux = true ;
        if ~quiet
            fprintf(1,'Solution is unique; rank condition met.\n') ;
        end
    else
        ux = false ;
        if ~quiet
            fprintf(1,'Multiple solutions exist; rank condition is not met.\n') ;
        end
    end
    
    exu = [ex,ux] ;
else
    exu = [] ;
end

% compute solution output form
%     k(t) = Tk + Tkk*k(t-1) + Tkw*w(t-1) + Tkv*v(t)
%     y(t) = Ty + Tyk*k(t) + Tyw*w(t)
Tk = gS ;
Tkk = Nss\Jss ;
Tkw = gSW ;
Ty = R*[zeros(ns,1); gU] ;
Tyk = R*[eye(ns); zeros(nu,ns)] ;
Tyw = R*[zeros(ns,nw); gUW] ;

[UU,SI,VV] = svd(Atilde) ; % UU*SI*VV' = Atilde ;
Tkv = VV*pinv(SI)*UU'*( C1*gamV - A1*Ru*gUW*phiV ) ;
end

function H = findvss(A,B,C)
% uses tricks with vec() and kron() to solve matrix equations of the form
%     H - A*H*C = B
% for the matrix H, given conformable matrices A, B, and C.
na = size(A,1) ;
nc = size(C,1) ;

H = reshape( (eye(na*nc)-kron(C',A))\vec(B), na, nc ) ;
end

function [L,q,tol] = span(A,B)
% Checks whether span(B) <= span(A) for two matrices A and B of dimensions
% m-by-na and m-by-nb. Follows the algorithm from Sims (2002)
% recapped in JKL.

tol = sqrt(eps) ;

U = opsvd(A) ;
T = opsvd(B) ;

n = size(U,1) ;

X = (eye(n)-U*U')*T ;

if norm(X)<tol
    L = true ;
    q = 0 ;
else
    L = false ;
    q = norm(X)-tol ;
end
end

function varargout = opsvd(A)
% [S,V,D] = opsvd(A)
% Outer product form of the singular value decomposition.

[S,V,D] = svd(A) ;

k = sum(sum(abs(V)>eps)) ;

varargout = cell(nargout,1) ;
varargout{1} = S(:,1:k) ;
if nargout>1
    varargout{2} = V(1:k,1:k) ;
    if nargout>2
        varargout{3} = D(1:k,:) ;
    end
end

end








